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We use computer simulation results for a dense Lennard- Jones fluid for a range of temperatures 
to test the accuracy of various binary collision approximations for the memory function for density 
fluctuations in liquids. The approximations tested include the moderate density approximation of 
the generalized Boltzmann-Enskog memory function (MGBE) of Mazenko and Yip, the binary colli- 
sion approximation (BCA) and the short time approximation (STA) of Ranganathan and Andersen, 
and various other approximations derived by us using diagrammatic methods. The tests are of two 
types. The first is a comparison of the correlation functions predicted by each approximate memory 
function with the simulation results, especially for the self longitudinal current correlation function 
(SLCC). The second is a direct comparison of each approximate memory function with a memory 
function numerically extracted from the correlation function data. The MGBE memory function is 
accurate at short times but decays to zero too slowly and gives a poor description of the correlation 
function at intermediate times. The BCA is exact at zero time, but it predicts a correlation func- 
tion that diverges at long times. The STA gives a reasonable description of the SLCC but does not 
predict the correct temperature dependence of the negative dip in the function that is associated 
with caging at low temperatures. None of the other binary collision approximations is a systematic 
improvement upon the STA. The extracted memory functions have a rapidly decaying short time 
part, much like the STA, and a much smaller, more slowly decaying part of the type predicted by 
mode coupling theory. Theories that use mode coupling commonly include a binary collision term 
in the memory function but do not discuss in detail the nature of that term. It is clear from the 
present work that the short time part of the memory function has behavior associated with brief 
binary repulsive collisions, such as those described by the STA. Collisions that include attractive as 
well as repulsive interactions, such as those of the MGBE, have a much longer duration, and theories 
that include them have memory functions that decay to zero much too slowly to provide a good first 
approximation for the correlation function. This leads us to speculate that the memory function 
for density fluctuations can be usefully regarded as a sum of at least three parts: a contribution 
from repulsive binary collisions (the STA or something similar to it), another short time part that is 
related to all the other interactions (but whose nature is not understood), and a longer time slowly 
decaying part that describes caging (of the type predicted by mode coupling theory). 

PACS numbers: 05.20.Dd, 05.20.Jj, 61.20.-p 



I. INTRODUCTION 

The concept of uncorrelated binary collisions^ has 
played an important role in the development of the ki- 
netic theory of fluids^ It is a central idea in the Boltz- 
mann kinetic equation for dilute gases, the Enskog equa- 
tion for high density liquids, and generalizations of En- 
skog theory that take into account the softness of the 
repulsive potential4^£&£ The nature of the binary col- 
lisions in the Boltzmann and Enskog theories are very 
different, involving in one case the entire interparticle 
potential, and in the other case only the repulsive part 
of the potential, which is idealized as a hard sphere po- 
tential. 

For kinetic theories of dense liquids that do take into 
account the effect of both the attractive and repulsive 
parts of the interatomic potential on the dynamics, the 
different range, strength, and effects of the two parts of 
the potential has motivated the construction of theories 
that treat these two parts very differently. Examples 
include the Rice-Allnatt theory^^ and the Karkheck- 



Stell theory^^^^^^? 

The equilibrium theory of the structure and thermo- 
dynamics of liquids shows similar developments. At 
low density, both the radial distribution function and 
the thermodynamic properties are represented by rel- 
atively simple expressions that contain the entire po- 
tential. At higher density, various theories have in- 
voked the primacy of the repulsive forces for determining 
the structure i£J&2LZL2LZL21 Theories that take into ac- 
count both types of forces at high density often calculate 
the effect of attractions using very different ideas from 
those used for repulsive forces i 18 i 21 ' 22 i 25 ' 26 

Theories of equilibrium structure and thermodynamics 
have been derived and formulated in a variety of ways, 
including perturbation theory, integral equations, and 
cluster expansionsj 27 ' 28 The Mayer cluster theory and its 
extension s 28 ' 29 ! 30 ' 31 provide a unifying theoretical frame- 
work for deriving and analyzing many of these theories 
(and many others) and relating them to one another. 
The Mayer theory represents the static correlation func- 
tions and free energy of a fluid in terms of an infinite 
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set of diagrams. Approximations are typically devel- 
oped by identifying certain subsets of diagrams, sum- 
ming them, and discarding the others. A standard tactic 
in cluster theory is to decompose the objects in the the- 
ory that represent interactions into various pieces and 
include the various pieces in different ways in developing 
approximations! 23 ' 24 ! 30 ' 31 

A focus on interparticle forces and their effects on dy- 
namics is quite different from the perspective provided 
by the fully renormalized kinetic theory of Mazenko and 
co-workers . 32 ' 33 In the latter approach, the the- 
ory of the correlation functions for density fluctuations 
in equilibrium liquids is formulated in such a way that 
the interparticle potential does not appear; instead, only 
the equilibrium static correlation functions, such as the 
pair correlation function and its generalization to more 
than two particles, appear in the final results. In effect, 
the dynamics is expressed in terms of the potentials of 
mean force rather than the 'bare' potential. 

One of the results of this theory was a generalized 
Boltzmann-Enskog memory function (denoted GBE) for 
fluids and an approximation to this function that they 
called the 'moderate density approximation' to the gen- 
eralized Boltzmann-Enskog memory function (denoted 
MGBE) . The latter can be interpreted physically as tak- 
ing into account uncorrelated binary collisions among the 
particles, with the potential that determines the dynam- 
ics being the potential of mean force rather than the bare 
potential. 

In recent work, Young and Anderse n 37 ' 38 used molec- 
ular dynamics simulations to study the behavior of sev- 
eral pairs of atomic liquids. Each pair consisted of two 
systems at the same density but different temperatures, 
with very similar pair correlation functions, and hence 
very similar potentials of mean force, despite their very 
different interparticle potentials. It was found that some 
(but not all) of the features of the dynamics of such a pair 
of liquids were very similar, suggesting that the potential 
of mean force is in some sense a more important determi- 
nant of dynamics than is the bare potential itself. This 
is consistent with the basic ideas of renormalized kinetic 
theory. 

A formally exact diagrammatic formulation of the ki- 
netic theory of density correlation functions was devel- 
oped by one of us i 39 ' 40 ' 41 It expresses the dynamic corre- 
lation functions of the particle density in terms of an infi- 
nite series of diagrams in a way analogous to the infinite 
series in the Mayer cluster theory for equilibrium static 
correlation functions. The vertices of the graphs, which 
represent interactions among particles, are expressed in 
terms of static correlation functions, and the theory is 
fully renormalized in the sense of Mazenko. This theory 
can provide a unifying theoretical framework for deriving 
and analyzing kinetic theories of density fluctuations in 
liquids. 

Using this graphical theory, Ranganathan and 
Anderse n 42 ' 43 considered the case of a fluid whose in- 
teratomic potential contained a continuous but very re- 



pulsive short ranged part, as well as longer ranged at- 
tractions. Using the inverse of the strength of the short 
ranged force as a small parameter, they determined the 
diagrams that are most important for the short time be- 
havior of the memory functions for the dynamic correla- 
tion function for density fluctuations. This analysis in- 
volved a procedure much like the one discussed above for 
the Mayer theory, namely a decomposition of the fully 
renormalized quantities that appear in the diagrams of 
the theory into various contributions with different char- 
acteristics in the limit as the small parameter goes to 
zero. The most important diagrams (i.e. the most di- 
vergent diagrams) were identified. Summing this infinite 
series of diagrams exactly led to an approximation for the 
memory function that is called the short time approxima- 
tion (STA). The structure of this approximation indicates 
that it takes into account uncorrelated binary collisions 
involving the repulsive part of the potential, with a colli- 
sion frequency that is determined by the structure of the 
fluid (and hence is influenced by the attractive forces). 

Using this short time approximation and neglecting the 
longer time contributions to the memory function gives 
a kinetic theory that can be regarded as a generalization 
of the Enskog theory to dense fluids with continuous re- 
pulsive forces as well as attractive forces. The Enskog 
kinetic equation is local in time, with no memory of the 
past history. As a result, the memory function in the 
Enskog theory is proportional to a Dirac delta function 
in time. For fluids with continuous repulsive potentials, 
there is a part of the memory function that is very large at 
very short times^ and that approaches a delta function 
in time as the repulsive part of the potential approaches 
a hard sphere potential. The STA includes all diagram- 
matic contributions to that large, short-lived part of the 
memory function. This is the sense in which the STA is 
a generalization of the Enskog theory. 

The predictions of this theory have been tested by com- 
parison with computer simulation studies of the Lennard- 
Jones fluid^ 3 - The accuracy of the theory is different for 
different correlation functions but is surprisingly good 
at high temperatures for the viscosity and diffusion con- 
stant. The theory fails to describe the changes in corre- 
lation functions that take place as a dense high temper- 
ature fluid is cooled to the triple point temperature. In 
particular, the increasingly negative dip in the velocity 
autocorrelation function and the self-longitudinal current 
are not well described, suggesting that the dynamics in 
the vicinity of the triple point requires much more than 
uncorrelated repulsive binary collisions for its descrip- 
tion. 

Ranganathan and Andersen also identified a more gen- 
eral set of diagrams whose sum leads to an approxima- 
tion, called the binary collision approximation (BCA), 
that describes the effect of uncorrelated binary collisions 
that are determined by the potential of mean force. (The 
relationship of this to the moderate density approxima- 
tion of the generalized Boltzmann-Enskog memory func- 
tion will be discussed in more detail below.) 
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In this paper, we discuss a variety of ways in which 
binary collisions might be defined for use in approxi- 
mations based on uncorrelated binary collisions. We 
test several of them by comparing their predictions with 
the results of molecular dynamics simulations of a dense 
Lennard-Jones liquid at various temperatures, from the 
triple point temperature to four times that value. There 
were several motivations for this work. 

1. There are many different possibilities for specifying 
the nature of the binary collisions in dense fluids. Some of 
them are mentioned above. Others are suggested by the 
possibility of breaking the potential of mean force, which 
appears in fully renormalizcd kinetic theory, into vari- 
ous contributions and describing the effect of the various 
contributions in different ways, in analogy to the way in 
which the bare potential is decomposed in both kinetic 
theory and the equilibrium theory of fluids. Moreover, 
the diagrammatic kinetic theory suggests some straight- 
forward extensions of binary collision approximations (as 
is discussed below). 

2. It would be worthwhile to test the various approx- 
imations by direct comparison of their predictions with 
computer simulation studies in order to determine which 
of them provides, in some sense, the best starting point 
for further development of kinetic theory. 

3. The mode coupling theory of Gotze and cowork- 
ers, and the kinetic theory that provides the basis for 
it , 45 i 46 ' 47 focusses on the long time behavior of the mem- 
ory function for density correlations in a fluid, but it as- 
sumes that there is a short time part that represents brief 
and presumably binary collisions between the atoms. 
However, the theory does not provide a microscopic ex- 
pression for this short time part, and in practice its con- 
tribution is described by a simple empirical function of 
time obtained by fitting to simulation or experimental 
data. It would be worthwhile to understand the nature 
of the true short time behavior of the memory function 
as a starting point for understanding the approximations 
made in mode coupling theory and attempting to go be- 
yond those approximations. 

Our discussion will rely heavily on the diagrammatic 
kinetic theor y 39 i 40 i 41 ' 42 as a unifying framework for stat- 
ing and comparing the various approximations as well 
as for deriving some of them. We restrict attention 
to dense atomic liquids (rather than gasses) with short 
range forces (rather than Coulombic forces). We also re- 
strict attention to theories that calculate and use a mem- 
ory function to calculate the equilibrium time correlation 
functions of the density in phase space. Finally, the cor- 
relation functions we focus on in this paper are self (or 
incoherent) correlation functions. 



II. BINARY COLLISION APPROXIMATIONS 
AND THE DIAGRAMMATIC THEORY 

A. Self correlation function 

We are concerned with the kinetic theory of the self 
correlation function of density fluctuations in a dense 
atomic fluid. The incoherent scattering function and self- 
diffusion are experimental observables associated with 
this function. (For additional studies that deal with other 
correlation functions, see Ref. l48l) 

The system of interest is a single component classical 
atomic fluid at equilibrium. Its Hamiltonian is 

N N 
i=l i<j=l 

The self part of the two-point time correlation function 
for density fluctuations in single particle phase space is 

^(ii / ,*) = ifew,t)*/ 4 (i',o)\ 

where p = N/V is the number density, 

f i (l,t)=6(R 1 -v i (t))6(P 1 -P i (t)) 

is the density in single particle phase space for particle i, 
and 

Sfi(l,t) = fi(l;t) - (fi(l;t)) 

is a density fluctuation. We use numbers like 1, 2, etc. 
to represent points in single particle phase space, e.g., 
1 — RiPi, 2' = R,2'P2'- Angular brackets denote equi- 
librium ensemble averages. 

The kinetic equation for this correlation function is 

= j dn! J d3M^(13;t-t')C^(3l',t') (1) 

where the self part of the memory function is the 

sum of two terms 

M (s) (ll',t) = M f (11', t) + M c w (ll',t) 

The flow term is 

M / (ll',t) = -^i-V R 5(ll / )^) 

TO 

We use the notation 6(11') = 5(R 1 -R V )6(P 1 -P V ) and 
Vr5(11') is the gradient of 6(11') with respect to the first 

(s) 

position argument. We make approximations for Mc , 
the collisional part of the self memory function, which we 
shall simply refer to as the memory function. 
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The initial condition for eq. (JTJ) is 

cW(ll',0) = pM ro (l)*(ll') 

where M m is the Maxwell-Boltzmann distribution of mo- 
menta 

M m (l) = (27rTOfc B T)~ 3/2 exp(-|P 1 | 2 /2mfc B T) 

and .g(ll') = ,g(|Ri — R-i'|) is the pair correlation func- 
tion of the fluid. 



B. Diagrammatic series 

The diagrammatic kinetic theory provides several 
closely related expressions for the correlation function 
C^ 8 ' and the memory function as the sum of the 

values of an infinite series of diagrams of a certain topo- 
logical structure i 39 i 40 ' 42 The diagrams contain propaga- 
tors that describe the time evolution of density fluctu- 
ations and vertices that represent renormalized interac- 
tions among those fluctuations. We shall focus on the 
two formulations of these series that are most useful for 
the present discussion. 

The first form of diagrammatic series 40 contains \^ 
propagators, which describe density fluctuations at sets 
of one or more points in the fluid, and Q vertices. 

The second form of the diagrammatic series 4 ^ contains 
-^(/p) propagators, each of which describes the propaga- 
tion of a single density fluctuation associated with a sin- 
gle particle that moves according to free particle motion, 
and Q( c,p ) vertices, which describe the localized interac- 
tions among such density fluctuations. 

The two series are each formally exact. The second is 
derived from first. The first is useful for proving some 
formal properties of the series, such as symmetry prop- 
erties. In practice, the second is more useful for deriving 
practical approximations. In the second series it is possi- 
ble to identify density fluctuations in different parts of a 
diagram and interpret them as being caused by the pres- 
ence (or absence) of the same particle. The vertices in 
the theory are calculated from static correlation functions 
of the fluid, and so it is possible to identify how many 
distinct correlated particles are associated with a vertex. 
The identification of the number of particles associated 
with a vertex is in some cases ambiguous, for reasons 
that we will discuss, but it is still possible to identify 
terms that correspond to, for example, the interaction of 
two particles, and hence it is possible to discuss whether 
diagrams describe binary collisions. 

Additional diagrammatic series can be derived from 
the second series by expressing the vertices as sums of 
various contributions. For example, in vertices that con- 
tain the potential of mean force, the latter can be de- 
composed into a short ranged repulsive part and a longer 
ranged oscillatory part, leading to the original vertex be- 
ing expressed as a sum of two new vertices. The resulting 



diagrammatic series can then be used to develop approx- 
imations in which the effects of the two parts of the po- 
tential are taken into account in different ways. (When 
an approximation is made to a diagrammatic series for a 
correlation function, by neglecting some diagrams and/or 
approximating certain vertices, the resulting approxima- 
tion may or may not be consistent with known symmetry 
properties of the correlation function, such as stationar- 
ity and time reversal symmetry. See Appendix [X] for a 
discussion.) 



C. Diagrammatic characterization of various 
approximations 

a. Generalized Boltzmann-Enskog memory function. 
This memory function, denoted GBE, is generated in 
Mazenko's theory by neglecting a function T, which rep- 
resents renormalized dynamical interactions, in the equa- 
tion of motion for a function in terms of which the mem- 
ory function of interest is expressed^ In terms of the 
first diagrammatic series, it represents the sum of all di- 
agrams in the series for the memory function in which 
there are at most two density fluctuations propagating 
at any time. 

This approximation satisfies important symmetry 
properties exactly (see Appendix [A]) . Moreover, it gives 
exactly the correct value of the memory function at t = 0. 
However, it is not clear how it might be calculated in any 
practical fashion, and to our knowledge numerical calcu- 
lations of the GBE memory function have never been 
performed. The physical meaning of the approximation 
is hard to specify but the mathematical meaning in the 
context of the diagrammatic theory is clear. It contains 
contributions to the memory function that describe two 
density fluctuations that interact with each other and 
propagate. The identities of the two particles associated 
with the two fluctuations can be different at different 
times. All interactions and propagations of the two fluc- 
tuations are treated exactly. The approximation that is 
made is to neglect the possibility of three or more simul- 
taneous fluctuations. 

b. Binary collision approximation. This approxima- 
tion, denoted BCA, was defined in Ref . l42l in terms of the 
second diagrammatic series. It is the sum of all diagrams 
in the series for the GBE memory function in which the 
same pair of particles is associated with the fluctuations 
at all times. 

This approximation does not satisfy important sym- 
metry properties exactly (see Appendix [A)) . However, it 
does give exactly the correct value for at t = 0. A 

practical algorithm exists to evaluate it, as is discussed 
below. 

The physical meaning of the approximation is that it 
contains contributions to the memory function that de- 
scribe two density fluctuations, associated with two spe- 
cific particles, that interact with each other and prop- 
agate, without having three fluctuations present at any 
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time. All interactions and propagations of the two fluc- 
tuations are treated exactly. The approximation here is 
not only to neglect the possibility of three or more simul- 
taneous fluctuations (as in the GBE) but also to neglect 
diagrams in which the two fluctuations at one time repre- 
sent a different pair of particles than the two fluctuations 
at other times. This characterization is the reason why 
it can be regarded as a binary collision approximation, 
i.e. one that involves only two particles. 

c. The moderate density approximation to the gener- 
alized Boltzmann-Enskog memory function. This mem- 
ory function, denoted MGBE, was defined by Mazenko as 
an approximation to the GBE. Various quantities in the 
GBE arc defined in terms of static correlation functions 
of the fluid. To construct the MGBE, the contributions 
to these quantities that contain three and four particle 
correlation functions are deleted^ 

The deletion of the three and four particle correlation 
functions to define the MGBE is closely related to, but 
not equivalent to, what is done in defining the BCA to 
include diagrams in which fluctuations representing at 
most two distinct particles appear. (See Appendix [Bl) 
The MGBE is almost equivalent to the BCA, except for 
one factor in every diagram. From a diagrammatic point 
of view, the MGBE can be regarded as an approximation 
to the BCA in which the vertex Q12 that appears in each 
diagram is replaced by an approximation that contains 
the potential of mean force. (A detailed expression for 
the MGBE is given below, together with a comparison to 
the BCA.) 

This approximation satisfies important symmetry 
properties exactly (see Appendix [A"]l . It does not give ex- 
actly the correct value of the memory function at t = 0. 
A practical algorithm exists to evaluate it, as is discussed 
below. 

The physical meaning of the MGBE is equivalent to 
that of the BCA. The only difference is that the addi- 
tional approximation made to the Q12 leads to a binary 
collision approximation that satisfies the symmetry con- 
dition. 

d. Short time approximation. This will be referred 
to as the STA. This was derived^ from the second di- 
agrammatic series using the procedure discussed above, 
namely expansion of the diagrams in powers of a small 
parameter representing the inverse of the strength of the 
repulsive part of the potential and retaining those terms 
that at short times are most divergent. The diagrams 
that appear are those of the BCA, as defined above, but 
with vertices that contain the repulsive part of the inter- 
atomic potential, rather than the potential or the poten- 
tial of mean force. (See below for a more precise state- 
ment of the value of the STA.) 

The STA satisfies the symmetry property approxi- 
mately, but not exactly. (In the limit of infinitely re- 
pulsive forces, it satisfies the symmetry property exactly. 
See Appendix [A]) It does not give exactly the correct 
value of the memory function at t — 0. A practical algo- 
rithm exists to evaluate it. 




FIG. 1: The first two diagrams are members of the diagram- 
matic series for the binary collision memory functions. (The 
remaining diagrams in the binary collision memory functions 
have additional Q^ 2 ' vertices on the two horizontal parti- 
cle paths.) The third and fourth diagrams are examples of 
diagrams that are not in a binary collision memory function. 
One or both of the two particles is interacting separately with 
a third particle in an STA approximation. The parts of the 
diagrams that have these STA interactions are indicated with 
dashed lines and dashed circles. (See the text and Fig. [2] for 
additional discussion.) 



e. Other binary collision approximations. Below we 
shall discuss several approximations whose diagrammatic 
meaning is straightforward. Like the STA, they can be 
defined graphically by starting with the BCA, separating 
each of the vertices into two terms, and retaining only the 
first term. The physical meaning of such an approxima- 
tion depends on the nature of the separation. The reason 
for investigating such approximations is the desire to de- 
velop better binary collision approximations as a starting 
point for a more general kinetic theory. 

/. Beyond the binary collision approximation. A 
major restrictive assumption of all the approximations 
mentioned above is that the possibility of having three 
or more simultaneous density fluctuations is ignored. The 
pair of particles that are associated with the two density 
fluctuations interact with each other, but the presence of 
other particles is taken into account only in some aver- 
age sense, if at all, due to the appearance of static pair 
correlation functions of the fluid in the various vertices. 

It is relatively simple, in the context of a graphical 
theory, to construct approximations that relax these as- 
sumptions in a physically meaningful way. The first two 
diagrams in a binary collision approximation are shown 
as the first two diagrams in Fig. [TJ The generic BCA 
diagram consists of two particle paths^ that propagate 
from the Q^i^ ver t ex on the right to the vertex 
on the left. In between these vertices arc zero or more 

(c 2) 

Q 2 2 vertices. These diagrams describe the propagation 
of a pair of particles that interact with each other using 

(c 2) 

the potential that is contained within the Q22 function. 
To develop the K approximations, consider including di- 
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FIG. 2: Examples of diagrams in a binary collision approxi- 

(s) 

mation augmented with Mg^ A memory functions on the par- 
ticle paths to represent interactions with third particles of 
the type shown in Fig. [1] The first diagram here includes the 
third diagram of Fig. [1] The second diagram here includes 
the fourth diagram of Fig. [T] The generic member of this se- 
ries is obtained by taking the generic member of the series for 
a binary collision approximation (see the first two diagrams 
in Fig. [T] for examples) and inserting an arbitrary number of 
Msta memory functions on each of the lines. (See the text 
and Fig. [T]for additional discussion.) 



agrams in the exact graphical series that are similar to 
these and that are constructed by inserting STA mem- 
ory functions into the propagators on the particle paths 
of the original diagrams of the BCA. The diagrams that 
result are those of the original BCA diagrams with any 
number of STA memory functions inserted. See Fig. [51 

When the diagrams with STA memory functions are 
evaluated, the two points have different times associated 
with them and the STA memory function is a function 
of the time difference. Since the collisions are brief, we 
approximate the effect of these vertices using a BGK 
model. 5 - The BGK model memory function is 

M K (ll';t) = vm K (ll')5(t) 

where m K (ll') = -£(Pi - Pi) + M m (Pi) and v is a 
collision frequency. The BGK model is an approximate 
memory function that is instantaneous in time. It ap- 
proximates a collision as having the effect of replacing the 
momentum of the particle by a new random momentum 
drawn from a Maxwell-Boltzmann distribution. Because 
of its simplicity, its effect can easily be incorporated into 
the algorithm for evaluating binary collision approxima- 
tions, as is discussed below. 

This procedure allows any binary collision approxima- 
tion to be converted into one that also includes brief un- 
corrected binary collisions of each of the two particles 
with third particles or fluctuations in the surrounding 
fluid. Graphically this is done by inserting Mjc vertices, 
as approximations for STA memory functions, into the 
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propagators that by themselves describe free parti- 
cle propagation. When this procedure is applied to any 
approximation, the resulting approximation will be given 
a name with "K" as a prehx. Thus, for example, the 
MGBE approximation, when modified in this way, will 
be called the K-MGBE approximation. In constructing 



such an approximation, the value of the parameter v must 
be chosen. This will be discussed below. 

We shall refer to this class of approximations as K ap- 
proximations. K approximations are mathematically and 
physically similar to a type of approximation suggested 
by Mazenko and Yip^ 



D. Analytic formula for the BCA and related 
appr oximat ions 

The various approximations mentioned above, with the 
exception of the GBE approximation, can be expressed 
in terms of very similar analytic formulas for M^ s \ As 
a result, all of them can be evaluated with a similar al- 
gorithm. Here we will state the formulas. The algorithm 
will be discussed later. 

For the BCA, the memory function is 

M ( gl A (lt,l't')= J d2d3did5 (2) 
x Q^ 1) (l;23)x S cA(23i;45t')Q^ 1) (45;l') 

where the Q vertices represent interaction terms and the 
Xbca term is a response function. The detailed equa- 
tions for these functions are 

Q^\l; 1'2') = V R T4 2 (12') ■ VpJ(ll') (3) 
Qi?\l2;l>) (4) 

= ^^p(2)M m (%(12)V R 7 21 (12) • VpJ(II') 

XBCA(l2t; 1'2't') (5) 
- e -iCW(t-t')xfn> 



6{ll')6(22')e(t - 1') 



The Liouville operator C is 
£(12) 



£(12) = L (l) + Lo(2)+Li(12) 
Pi 



m(l) 



• V 



Ri 



(6) 
(7) 



Li(12) = i V R F 22 (12).(V Pl -Vp 2 ) (8) 

These formulas contain three functions V12, V21 and V22 
that are two particle potentials. For the BCA, they are 

V 12 (12) = u(12); V 21 (12) = V 22 (12) = v MF {12) (9) 

where u{12) is the interparticle potential and vmp(12) = 
— fc^Tln (?(12) is the potential of mean force, and g(12) 
is the pair correlation function. 

The BCA is exact in the low density limit and corre- 
sponds to the linearized Boltzmann equation in that limit 
on appropriate length and time scales. The xbca de- 
scribes two particle scattering with an effective interpar- 
ticle potential that is the potential of mean force, which 
also appears in Q 2 \. At low density, the potential of 
mean force approaches the bare potential. The BCA de- 
scribes binary collisions (on microscopic length and time 
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Approximation 


V12 


V21 


V22 




STA 


Ur(r) 


u r (r) 


u r (r) 


-332 


BCA 


u(r) 


vmf (r) 


VMF(r) 


-281 


K-BCA 


u(r) 


VMF(r) 


VMF(r) 


-281 


MGBE 


vmf{t) 


vmf (r) 


vmf(t) 


-289 




VMF\f) 


vmf{t) 


vmf{t) 


-zoy 


RPMF 


V r Mp{r) 


VrMF(r) 


VrMF{r) 


-215 


K-RPMF 


v,-mf(j) 


VrMF(r) 




-215 


hybrid 
MGBE/STA 


vmf{v) 


vmf{t) 


u r (r) 


-289 



TABLE I: A table of all approximations studied, stating the 
necessary potential functions that define a binary collision 
approximation for the memory function and the value of the 
kk matrix element of the memory function at t = 0. The 
BCA value for M^p(k, 0) is exact. Note that u(r) is the bare 
potential, u r (r) is the repulsive part of the bare potential, 
vmf(j) is the potential of mean force, and v r MF(r) is the 
repulsive part of the potential of mean force. Details about 
approximations not discussed in this text can be found in 
Ref . IH 



III. CALCULATIONS 
A. Basis function expansion of the kinetic equation 

Using the expressions for the Q vertices and xbca re- 
sponse function, the expression in eq. ([2]) leads to 

<a(«; 1'*') = / d2d3dAd5^^- )P (5)M m (5)g(45) 

x (V R tq 2 (13) • V P <5(12)) ■ (V R T/ 21 (45) • V P <5(41')) 
x e l£(23)(t -*' ) (5(24)5(35)e(i - t') 

It is convenient to take the Fourier transform of C^ s ' 
and with regard to the difference in their position 

arguments and to represent their momentum dependence 
in terms of a set of Hcrmitc polynomial basis functions. 
The details of the Hermite polynomial expansion are dis- 
cussed in Ref. |48|. From eq. ©, we can obtain the fol- 
lowing expression for the Fourier transformed matrix el- 
ements of 



scales) in much the same way as the Boltzmann equa- 
tion, with some many body effects accounted for in the 
effective potential. 

The MGBE is of the same form as eqs. ©h®, with 
© replaced by 

V 12 (12) = V 21 (12) = V 22 (12) = v M f{12) (10) 

The STA is of the same form as eqs. ©-©, with © 
replaced by 

V 12 {12) = V 21 {12) = V 22 {12) = u r (\2) 

where u r {\2) is the repulsive part of the potential it(12). 

Any other binary collision approximation, as discussed 
above, leads to results of the form of eqs. with 
the potentials Vi 2 , V 2 \, and V 22 being the functions used 
in the approximation. See Table U for a list of various 
approximations we have considered. 

Any K approximation, as discussed above, leads to 
results of the form of eqs. ([2])-© with eq. © replaced 
by 

C K {12) (11) 
= L (l) + Lo(2) + Li(12) + L BGK {1) + L BGK {2) 

The symmetry properties of the various approximations 
are discussed in Appendix [A"l 



M^bca^ t) = -J dp i dR 2 dP 2 <?(|Ri - Ra|) 

x M m (Pi)pM m (P 2 )(VVQ 21 (|Ri -R 2 |) • VMPi)) 
x e^ k - R »(VVQ 12 (|R a -R h |)-V^(P a )) (12) 

The jl and v are labels for the Hermite polynomial basis 
functions. Each is an ordered triplet of nonnegative inte- 
gers. (In the following, we shall use '0' as an abbreviation 
for the triplet '000' and k as an abbreviation for '001'. 
Note that M ( -f) = if either u = or v = 0. Mi 8 ) is 

f"^ " kk 

the single most important matrix element of the memory 
function. In the limit of small wave vector and/or short 
time, it is the only matrix element that contributes to 
the correlation function, and in general it is the matrix 
element that most strongly affects the time dependence 
of the correlation function.) In eq. (fT2")l . Ri is defined 
to be the origin and the phase points R1P1R2P2 are 
the initial conditions for two particles that evolve for- 
ward in time according to Hamilton's equations of mo- 
tion with V 22 as the interatomic potential to the final 
phase points R a P a Rf,P;, at time t. The matrix elements 
of the memory function are calculated using the two par- 
ticle trajectory calculation method of Ranganathan and 
Andersen42 See Appendix [C] for details and for a discus- 
sion of the error analysis. 

The correlation functions of interest in this work are 
the self correlation functions: the incoherent intermedi- 
ate scattering function (IISF) 

and the self longitudinal current correlation function 
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(SLCC) 

4(M) = (*»,*) 

\ I P»z W feW e -ifc a (rj,(t)-ri 8 (0))\ 

AT mm I 

In these equations it is to be understood that the wave 
vector k is in the z-direction. The IISF and SLCC func- 
tions are related in the following way 

4(k,t) = -^J^F s (k,*) 

In the graphs below, the SLCC is normalized so that its 
value at unity is t = 0. We focus on the SLCC because we 
have found its qualitative features to be sensitive to errors 
in binary collision approximations. For a discussion of 
additional correlation functions, see Ref. |48|. 

For some of the discussion we will need two related 
functions. The first is the total correlation function 

C(ll',t) = (Sf(l,t)5f(l',0)) 

where 

N 

/(l,t) = 5^5(R 1 -r i (t))«(P 1 -p i (t)) 

is the total particle density and 

Sf(l,t) = f(l,t)-{f(l,t)) 

is a density fluctuation. The second function is the 
Fourier transform of C(ll',i), the coherent intermediate 
scattering function (ISF) 

F(k,i) = <? W (M) = jf(jt e-* (riW - r ' (0)) ^ 



The IISF, which corresponds to C$ (k, t) , is converged 
with the 0(t 9 ) approximation. The SLCC, which corre- 
sponds to C^(k, t), is converged with the 0(t 7 ) approx- 
imation. We present only these solutions in the results 
that follow. Moreover, we concentrate on the results for 
the SLCC since none of the approximations are able to 
consistently reproduce all of the features of this function, 
especially the temperature dependence of the negative 
region for the function. Since the SLCC is proportional 
the second derivative of the IISF, the accuracy of an ap- 
proximation for the SLCC function is consistent with its 
accuracy for the IISF. 



C. Molecular dynamics calculations 

The correlation functions calculated from the binary 
collision approximations for the memory function will 
be compared to results generated from molecular dy- 
namics simulations of a 500-particle, one component 
Lennard-Jones liquid at equilibriuni i 37 i 38 i 55 The states 
of interest are reduced density p = 0.85 (approximately 
the triple point density) at reduced temperatures T — 
0.723, 1.554, 3.000 (the triple point temperature is ap- 
proximately 0.723). Results for temperatures between 
0.723 and 1.554 at p = 0.85, and high temperature results 
at p = 0.75 can be found in Ref. |48|. The usual reduced 
units for the Lennard-Jones fluid are usedi^ The corre- 
lation function results from molecular dynamics simula- 
tion and the binary collision approximations have small 
statistical error whose values are generally smaller than 
the differences between the simulation data and the the- 
oretical curves. The details of the molecular dynamics 
simulations can be found in Refs. [37| and l38l . 



B. Solution of the kinetic equations 

Use of the Hermite polynomial basis set for the mo- 
mentum arguments converts the kinetic equation for 
(k, t) into an infinite dimensional matrix equation 
for the matrix elements of & s >. We solve this infinite 
dimensional equation numerically using the method of 
Ranganathan and Andersen^ 3 - In this method, the equa- 
tions for successively larger finite sets of coupled equa- 
tions are solved. Each of these finite sets is chosen such 
that its solution agrees with the exact solution of the in- 
finite set up to order t n for small t, and successive finite 
sets have increasing values of n. When the successive so- 
lutions no longer change with increasing n, the last result 
is regarded as the solution to the infinite set. A detailed 
discussion of the method can be found in Ref. [54l. 



IV. RESULTS 
A. BCA 

Although the BCA approximation gives exactly the 
correct t = value of the memory function, this ap- 
proximation leads to correlation functions that oscillate 
rapidly and with sharply increasing magnitude as t — > oo. 
Although the shape of the BCA memory function looks 
reasonable, there is insufficient area under the curve for 
the solution to the kinetic equation to be well-behaved. 
We suspect that this is related to its failure to satisfy the 
symmetry property discussed in Appendix[XJ (Recall the 
form of this approximation, eqs. ©-{H]).) 
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B. STA, MGBE, and K-MGBE 

1. STA 

Overall, for the high densities studied, the STA ap- 
proximation does reasonably well at predicting correla- 
tion functions at all temperatures and wave vectors, with 
increasing quantitative accuracy with high temperature 
and large wave vector. The short time predictions of the 
STA are reasonably accurate, as is the overall longer time 
behavior. The STA approximation is able to reproduce 
qualitative trends of the various correlation functions 
with changes in wave vector. The trends with changes 
in temperature are also reproduced, except for the SLCC 
curve where the shifting and increasing magnitude of the 
negative dip is not well reproduced. In fact, the STA ap- 
proximation consistently fails to predict the temperature 
dependence of the negative dips in the SLCC curve. A 
complete discussion of the results of the STA approxima- 
tion can be found in Refs. and EH 

2. MGBE 

The MGBE approximation is more accurate at short 
times than the STA, as shown in Fig. [3J and is in es- 
sentially perfect agreement with the simulation data for 
times less than about 0.1. This statement holds for the 
self correlation functions at all temperatures and wave 
vectors studied. The short time part of the SLCC func- 
tion is largely determined by the value of the kk matrix 
element of the memory function at t = 0. The MGBE 
does not have exactly the correct zero time values of the 
memory function matrix elements but it is more accurate 
than the STA, and this accounts for its accuracy at short 
times. The MGBE result begins to become inaccurate at 
t w 0.1, rising significantly above the simulation data. 

For longer times, t > 0.1, the MGBE gives very inac- 
curate results. See Figs. 2] and [5] For large wave vectors 
at both high and low temperatures, the MGBE behaves 
reasonably, with a single minimum in the SLCC at ap- 
proximately the right time, but an overall quantitative 
disagreement with the data. For small wave vectors at 
all temperatures, it predicts oscillation in the SLCC for 
times of about 0.2-0.3, a positive maximum in the func- 
tion at t Ri 0.3, and a very slow decay of the correlation 
function to zero from above, in qualitative disagreement 
with the simulation results. The latter have a single min- 
imum, followed by an approach to zero from below on 
a shorter time scale than the MGBE result approaches 
zero. 

The origin of this unusual behavior, especially the os- 
cillations, is easily understood from the nature of the 
MGBE. The memory function in the MGBE is a time cor- 
relation function for derivatives of the potential of mean 
force in a two particle system for which the two particles 
interact according to the potential of mean force of the 
fluid and move according to Hamiltonian dynamics. See 
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FIG. 3: Normalized SLCC function. Although the MGBE 
and K-BCA approximations do poorly at representing the 
simulation result for longer times, for times t < 0.1 we see 
that these are the most quantitatively accurate approxima- 
tions. The statistical error bars for the MGBE correlation 
function are shown and are barely discernable on this scale. 
All of the theoretical correlation functions have similarly small 
statistical errors. 

eqs. ©, ([T0|). |12]) and Appendix [C] This idea, which 
follows from the explicit formula eq. (fT2|) , is the basis for 
the algorithm used to evaluate the matrix elements of the 
memory function. The potential of mean force at high 
density is oscillatory and extends to distances of several 
atomic diameters. Thus, collisions that are relevant for 
the MGBE are of long duration and the correlation func- 
tion that gives the memory function is nonzero for a very 
long time. 

Figure [6] shows the time dependence of the principal 
matrix element of the memory function for the MGBE 
approximation for a low temperature state. There is a 
pronounced oscillation that persists for times longer than 
the time scale for the true correlation function to decay. 
For comparison, the STA memory function decays to zero 
very quickly, within times of the order of 0.1. As we shall 
discuss below in the next section, the large positive fea- 
ture in the memory function and the subsequent oscilla- 
tion are not present in the actual memory function. 

These results show that the MGBE approximation, al- 
though quantitatively correct at short times, is qualita- 
tively incorrect at intermediate and long times. 



3. K-MGBE 

The collisions described above that are the basis of the 
MGBE approximation, namely Hamiltonian dynamics of 
a pair of particles whose interaction potential is equal to 
the potential of mean force of the liquid, seem inappro- 
priate for determining the dynamics of a dense liquid. 
The particles surrounding a pair of particles should do 
more to determine the dynamics of the particle pair than 
merely to renormalize the bare potential into the poten- 
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FIG. 4: Normalized SLCC function for the STA, MGBE and K-MGBE approximations for T = 0.723 at small and large wave 
vector. The STA and MGBE curves lie on top of each other at large wave vector. High temperature results are shown in Fig.[S] 
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FIG. 5: Normalized SLCC results for the STA, MGBE, and K-MGBE approximations for T = 1.554 at small and large wave 
vector. See Fig. U for low temperature results. 
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FIG. 6: Self part of the kk matrix element of the memory 
function. The STA approximation generates a non-oscillatory 
memory function while both MGBE and K-MGBE approxi- 
mations yield oscillatory memory functions. These oscilla- 
tions are the cause of the oscillations in the MGBE and K- 
MGBE results for the SLCC. 



tial of mean force. In Mazenko's theory, this would be 
taken into account by including the quantity called T, 
which represents renormalized dynamical interactions of 
the pair. In the cluster theory, this would be taken into 
account by including the coupling of the pair to fluctua- 
tions of the density of third particles. 

The K-MGBE approximation discussed above provides 
a tractable way of including these interactions, provided 
they are regarded as independent uncorrelated repulsive 
collisions of the two particles with third particles, approx- 
imated using the BGK model memory function. This 
approximation sums, in an approximate way, a class of 
diagrams in the exact diagrammatic theory. We assume 
that these collisions involve the short range repulsive part 
of the interatomic potential, so we pick the parameter v 
so that the BGK memory function has a strength that is 
related to that of the STA memory function. Since the 
most important matrix element of the memory function 
is the kk matrix element, we choose v so that this matrix 
element for the BGK model corresponds to the actual 
matrix element for the STA. More precisely, we choose v 
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such that 

/>oo 

"=/„ dt ' HisTA^ ( 13 ) 
«/ u 

While this is a crude approximation overall, it was de- 
vised to give some idea of the effects of third particles 
on the dynamics of the pair, in the hope of yielding an 
improved version of the MGBE. 

Figure [5] shows that the K-MGBE memory function 
has less pronounced oscillations and a more rapid decay 
to zero than the MGBE memory function. Figures [4] and 
[5] show the SLCC correlation functions given by the K- 
MGBE approximation. Introducing the BGK collisions 
extends the initial time interval over which the approx- 
imation is valid, for low and high temperatures and for 
small and large wave vector. The intermediate time be- 
havior of the K-MBGE is unsatisfactory, however, with 
a peak in the correlation function at times of about 0.3 
appearing under all conditions studied. The long time 
behavior of the K-MGBE is at best a slight improvement 
over that of the MGBE. 



C. Additional approximations 

In an attempt to find improved binary collision approx- 
imations, we investigated a number of other possibilities. 

a. RPMF approximation. The STA memory func- 
tion, which uses the repulsive part of the potential in 
each of the vertices, has the virtue of decaying to zero 
rapidly without oscillations at long time that produce 
unphysical behavior of the correlation function at long 
times. The rapid decay is a result of the fact that the 
two particle dynamics in the calculation of the memory 
function involves purely repulsive interactions (and hence 
brief collisions) and the range of the repulsive potential 
is so short. 

An approximation, which we denote RPMF, based on 
using the repulsive part of the potential of mean force 
in each vertex, was constructed and studied. While the 
memory function generated by this approximation has 
the feature of rapid decay, as discussed for the STA, its 
magnitude and its integrated area are smaller than that 
of the STA and it generates an SLCC that falls from 
its initial value significantly more slowly that the STA 
result. This basic flaw in the RPMF approximation is 
evident in the zero time value of the memory function 
(see Table HI). 

b. K-RPMF approximation. The K-RPMF approx- 
imation adds BGK collisions with third atoms to the 
RPMF approximation. This memory function decays 
more slowly than that of the RPMF, but it has the same 
initial value as the RPMF, and overall is not a significant 
improvement over the RPMF. 

c. A hybrid of the MGBE and STA approximations. 
We attempted to construct approximations that have ac- 
curate values of the memory function at zero time but 
that were based on relatively brief collisions. In one such 



approximation, we used the potential of mean force in the 
Q^vi^ vertex and the Q^i^ vertex (as in the MGBE ap- 
proximation) but the repulsive part of the potential in the 
Q22 vertex (as in the STA). We also considered the ap- 
proximation that supplemented this approximation with 
BGK collisions with third particles. This hybrid is not a 
significant improvement over the MGBE or STA approx- 
imations and had some artifacts for small wave vector 
that are inconsistent with the simulation results. 

d. Other approximations. A variety of other unsuc- 
cessful approximation were also studied. See Ref. [HI for 
details. 



V. EXTRACTION OF THE MEMORY 
FUNCTION 

Thus far, we have tested approximate kinetic theories 
by calculating their memory functions, solving the mem- 
ory function equation for the related correlation func- 
tions, and comparing the predicted correlation functions 
with computer simulation results. These comparisons 
give only an indirect test of the accuracy of the approxi- 
mate memory function approximations. 

In this section we consider a different set of questions. 
Given a correlation function calculated from simulation 
data, what is the memory function that corresponds to 
it? What are the properties of the 'exact' memory func- 
tion? How do the approximate memory functions com- 
pare with the 'exact' memory function? 

The numerical technique we use for extracting the 
memory function from the correlation function data em- 
ploys a Fourier series representation of the memory func- 
tion as a function of time. The Fourier coefficients in 
the series were adjusted to minimize the differences be- 
tween the SLCC calculated from the memory function 
and the SLCC obtained from computer simulation. For 
small wave vector, the kk matrix element of the mem- 
ory function plays the dominant role in determining the 
SLCC function^ The reason for this is that the cou- 
pling of the correlation function to other elements of the 
memory function involves positive powers of the wave 
vector. Thus we need to construct a Fourier representa- 
tion of only this matrix element of the memory function 
and solve for the Fourier coefficients of just one func- 
tion rather than many. The coefficients are found by a 
minimization of the sum of the squares of the deviations 
of the calculated correlation function and the simulated 
function. See Ref. |48| for details. 

The extracted (k, i) for small wave vector is shown 
in Fig. [7| for the lowest and highest temperatures stud- 
ied, where they are compared with the STA predictions 
for the same function. The extracted memory function 
at low temperature has none of the oscillations and high 
maxima that are predicted by the MGBE and K-MGBE 
approximations. Compare Fig. [5] The STA results are 
qualitatively similar to the extracted results, but there 
are significant differences. 1. The magnitude at zero time 
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FIG. 7: Extraction results (dashed lines) for the kk matrix 
element of the self memory function compa STA curves (solid 
lines) at p = 0.85 and |k| = 0.75. The top set of lines at t = 
corresponds to T = 0. The lower set at t = corresponds to 
T = 3.000. 



FIG. 9: The longer time part of the memory function pre- 
dicted by mode coupling theory is compared to the memory 
function numerically extracted from the simulation data for 
the SLCC at T = 0.723 and Ikl = 0.75. 




FIG. 8: The difference between the extracted and STA mem- 
ory functions, Af fcfciSTj4 (0.75, t) - M k k,ext{0-75,t), for three 
different temperatures. 



is larger for the STA than for the extracted memory func- 
tion. 2. The STA curves fall rapidly to zero, while the 
extracted curves decay to zero more slowly, from above 
at the higher temperature, and from below at lower tem- 
peratures. 

The differences between the extracted memory func- 
tion and the STA approximation is plotted for high, 
medium, and low temperatures in Fig. [8l At the highest 
temperature, the memory function decays to zero from 
above. This is presumably a manifestation of the hydro- 
dynamic vortex effect that is seen for hard spheres^ The 
long time negative part of the memory function at lower 
temperatures is presumably a manifestation of the cage 
effect, and is related to the long slow decay of the simu- 
lation correlation function at low temperature and small 
wave vector (see Fig. 0]). The very long time part of the 
memory function is shown in Fig. [5] 

The actual memory function at low temperatures 



clearly has a rapidly decaying part, falling to 2% of its 
initial value at t « 0.5, and a more slowly decaying part 
that falls to zero (within the statistical error) at about 
t = 1.5. None of the binary collision approximations we 
have studied is qualitatively consistent with the true long 
time behavior of the memory function. Moreover, none 
of them provides a quantitative description of the initial 
rapid decay. 



VI. MODE COUPLING DESCRIPTION OF THE 
LONG TIME PART OF THE MEMORY 
FUNCTION 



The mode coupling theory of Gotze^ can be used 
to calculate, among other things, a correlation func- 
tion denoted (j^ (t) that is proportional to (k, i) , 
the Fourier transform of a specific matrix element of 
C^ s > (11', t). The theory is based on the Mori-Zwanzig 
projection operator method,^ 9 " which relates this function 
to its second order memory function, denoted Mq(t). 
The assumptions of the theory lead to a mode coupling 
relationship that expresses the long time part of Mq S ^ (t) 
as a bilinear functional of 4> q s ^ (t) and 4> q (t) , where the lat- 
ter is proportional to Cjjjj(q, t). See eq. (3.34) of Ref. l58l. 

The kinetic theory presented here relates (k, t) 
to an infinite dimensional matrix of memory functions 
Mi s J (k, i). For small wave vector k, (ji s J (k, t) is coupled 

primarily to Mj- s ~ (k, t) , and in the limit as this wave vec- 
tor goes to zero, this is the only element that affects the 
value of C^g (k, t). Our numerical work confirms that 
k| = 0.75 is small enough for this to be the case. It 
is straightforward to show that when this is the case, 
the relationship between M- 8 - (k, t) and Ci s J (k, t) is pre- 
cisely of the form generated by the Mori-Zwanzig theory. 



13 



In other words, for small wave vector, the single ma- 
trix element (k, t) is the Mori-Zwanzig second order 

memory function of C^(k,t). 

When the mode coupling relationship mentioned above 
is expressed in terms of the present kinetic theory, the 
result is of the following form 

M« (k, t) = J dqA(k, q)C^(q, t)C$ (k - q, t) 

where A(k, q) is a nonnegative function. The mode 
coupling theory asserts that this holds for times longer 
than the duration of brief binary collisions, which in the 
present theory corresponds roughly to the time at which 
the STA memory function has decayed to zero. 
Applying this result to |k| = 0, we find 

(0,t) = Jdq A(0, q)C m (q, t)C$ (q, t) (14) 

Each of the two time dependent functions on the right 
depends on the magnitude, but not the direction, of its 
wave vector argument. The function A(0,q) is zero if 
q| = 0, small near |q| = 0, and highly peaked near 
q| = Qmax ~ 27r/<7. Assuming that the integral above is 
dominated by the range of wave vectors whose magnitude 
is near this value, we find 

M${0,t) « (constant) xC m (q max ,t)C^(q max ,t) (15) 
where the constant is positive. Since Ml s ) (k, t) is an even 

kk v ' 

differentiable function of the wave vector argument, it is 
very slowly varying near |k| = 0, and we will apply this 
result for |k| = 0.75, the smallest nonzero wave vector 
studied in this work. 

We can test this relationship by comparing the time 
dependence of the extracted M^(0.75,i) obtained from 
the simulation data with the time dependence of the right 
side of eq. (|15p . which contains functions that are rou- 
tinely obtained from simulation data. Figure [5] shows 
the extracted memory function My (0.75, t) for longer 
times (t > 0.1) for the lowest temperature studied. At 
that time, the STA memory function has decayed almost 
completely to zero. The extracted memory function has 
oscillations that are an artifact of the extraction proce- 
dure (note the vertical scale compared with the graph of 
Fig. [7]). Figure[9]also shows the right side of eq. lfl5|) . with 
the constant chosen so that the results matches the value 
of the extracted memory function for t w 1. With this 
choice of the constant, the time dependence of the right 
side accounts remarkably well for the time dependence of 
the memory function matrix element. 

Note that this test was only for the smallest wave vec- 
tor at the lowest temperature studied. For high temper- 
atures, eq. (fTS"|) is qualitatively incorrect (i.e., the pre- 
dicted sign of the memory function is incorrect). This is 
presumably because the behavior is dominated by a hy- 
drodynamic vortex effect^ rather than the cage effect 



described by mode coupling theory. For higher wave vec- 
tors at low temperatures, the test cannot be performed, 
since the simple relationship between My (k, t) and the 
second order memory function of the Mori-Zwanzig the- 
ory does not hold. Note, also, that mode coupling the- 
ory allows, in principle, the calculation of the constant of 
proportionality, but we have not done this. Nevertheless, 
the agreement in Fig. [5] is striking, suggesting that mode 
coupling theory accounts quantitatively for the time de- 
pendence of the memory function for long times (i.e., for 
times longer than the time of a brief binary collision) for 
the lowest temperature studied, which corresponds to the 
liquid state near the triple point. 

VII. SUMMARY AND CONCLUSIONS 

a. The accuracy of various approximate theories. 
We have investigated various binary collision approxima- 
tions for the memory function of the correlation function 
of density fluctuations in a dense Lennard- Jones liquid 
at equilibrium and compared their predictions with com- 
puter simulation results for the self correlation functions, 
especially the self longitudinal current correlation func- 
tion (SLCC). 

The BCA of Ranganathan and Andersen is exact at 
t = 0, but it leads to divergent results at long times, pre- 
sumably because of its failure to satisfy important sym- 
metry properties. 

The moderate density approximation of the general- 
ized Boltzmann-Enskog memory function (MGBE) of 
Mazenko and Yip2& is accurate at very short times. How- 
ever, its intermediate time behavior is very different from 
the correct results, and its predicted correlation functions 
display unphysical oscillations at intermediate times. 

The STA of Ranganathan and Anderse n 42 ' 43 is not 
exact at t = 0. It describes the initial drop of the 
SLCC to zero reasonably accurately, but it does not ac- 
curately describe the shape and depth of the minimum in 
that function and its temperature dependence. At high 
temperatures, it gives accurate values of the transport 
coefficients 4 s - 

The K-MGBE approximation, which introduces uncor- 
related brief collisions of the colliding pair with third par- 
ticles, is more accurate than the MGBE at short times 
but its intermediate time behavior is still very inaccurate. 

No other binary collision approximation discussed in 
this paper is a systematic improvement over the STA. 

b. The behavior of the extracted memory function. 
For small wave vector, the most important matrix ele- 
ment of the memory function can be extracted from the 
simulation results for the SLCC. 

At short time, the extracted memory function drops 
rapidly to values close to zero, in a way similar to the STA 
result, and then has a very small, long-lived, nonoscilla- 
tory tail that is positive at high temperatures and nega- 
tive at low temperatures. At low temperatures, the neg- 
ative tail is important for determining the shape of the 
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minimum in the SLCC. This tail is probably the result of 
physical effects that are not included in any binary col- 
lision approximation. The time dependence of the long 
time tail of the memory function is in striking agreement 
with what is predicted by mode coupling theory, but this 
is only a limited confirmation of that theory since the ex- 
tracted memory function can be obtained only at small 
wave vector. 

c. Conclusions and speculations. The negative fea- 
ture of the SLCC for the dense Lennard-Jones fluid in- 
dicates that the velocity of a particle is anticorrelated 
with its velocity at times that are more than about 0.1 
time units earlier. An important part of the basis for 
this feature can be described by binary collision theories 
that take into account the finite duration of very repul- 
sive collisions. For hard spheres, the short time feature 
of the memory function is a Dirac delta function, which 
by itself would generate monotonic decay of the SLCC 
from above. In the STA, the short time feature is spread 
over a range of times of the order of the duration of a 
repulsive collision, and the resulting correlation function 
has a negative dip. However, the increase of the depth 
of the negative feature as the temperature is lowered is 
caused not by this short time feature but rather by a 
small, negative, temperature dependent long time tail in 
the memory function. This is presumably a manifesta- 
tion of what is commonly regarded as caging of atoms at 
low temperatures and cannot be described by binary col- 
lision approximations. Instead, mode coupling theory or 
something else that goes beyond binary collision approx- 
imations is needed to describe this aspect of the negative 
feature. 

This conclusion is not at all surprising. Theories that 
use mode coupling commonly include a binary collision 
term in the memory function but do not discuss in de- 
tail the nature of that term. It is clear from the present 
work that the short time part of the memory function has 
behavior associated with brief binary repulsive collisions, 
such as those described by the STA. Collisions that in- 
clude attractive as well as repulsive interactions, such as 
those of the MGBE, have a much longer duration, and 
theories that include them have memory functions that 
decay to zero much too slowly to provide a good first ap- 
proximation for the correlation functions at intermediate 
times. 

We have not been able to derive any binary collision 
theory that describes brief binary repulsive collisions and 
that has exactly the correct value of the memory func- 
tion at t = 0. There is a straightforward expression for 
the zero time value of the memory function^ in terms of 
static correlation functions, and this expression involves 
the entire potential, the entire potential of mean force, 
and the pair correlation function. Any graphical approx- 
imation for the memory function that includes all this at 
t = would have more physics in it than binary repulsive 
collisions. 

This leads us to speculate that the memory function for 
density fluctuations can be usefully regarded as a sum of 



at least three parts: a contribution from repulsive binary 
collisions (the STA or something similar to it), another 
short time part that is related to all the other interactions 
(but whose nature is not understood), and a longer time 
slowly decaying part that describes caging (of the type 
predicted by mode coupling theory). 
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APPENDIX A: SYMMETRY PROPERTIES OF 
THE MEMORY FUNCTION 

The theory we are discussing describes a system with a 
time independent Hamiltonian at equilibrium. The exact 
correlation function for the system must satisfy symme- 
try properties related to stationarity and time reversal 
invariance. Stationarity implies 

C (s) (Ri , Pi ; R'i , Pi ; i) = C (s) (R'i , P'i ; Ri , Pi ; -t) 

. ( A1 ) 

A combination of time reversal symmetry and stationar- 
ity implies 

C« (Rx, P l5 Ri, Pi ; t) = (Ri , -Pi; Rr, -P l5 1) 

There are similar implications for other correlation func- 
tions and for the memory function. The importance of 
related symmetry conditions formulated in the Laplace 
transform domain have been emphasized by Mazenko. 36 

The exact diagrammatic series satisfies these relation- 
ships exactly. This follows from the fact that the se- 
ries are formally exact, but it can also be proven di- 
rectly from the topological structure of the diagrams 
and the symmetry properties of the vertices and bonds 
that appear in the diagrams. (The diagrammatic theory 
was derive d 39 i 40 ' 41 i 42 ' 43 only for positive time arguments. 
However, the results can easily be extended to negative 
time by analytic continuation in order to verify that they 
satisfy eq. (|ATjl .) 

Graphical approximations are constructed by retain- 
ing a subset of diagrams in the exact series and/or by 
making approximations for the vertices. Such graphical 
approximations may or may not satisfy the symmetry 
requirements, depending on the nature of the subset re- 
tained and details of the approximations made. 

We have not been able to formulate a useful set of 
necessary conditions that a general graphical approxi- 
mation must satisfy in order to be consistent with the 
symmetry properties. However, various sets of sufficient 
conditions can be formulated and proven. In particular, 
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two sets of sufficient conditions for a binary collision ap- 
proximation, of the class discussed above, to satisfy the 
symmetry property can be stated. The sufficient condi- 
tions reduce to conditions on the quantities that appear 
in the vertices. The details will be omitted here and we 
shall merely quote the results. 1. The MGBE approxima- 
tion for Af( s ) satisfies the symmetry conditions exactly. 
2. The BCA approximation for does not satisfy the 
known sufficient conditions for symmetry. It is likely that 
it strongly fails to satisfy the symmetry, and this is prob- 
ably the reason why it is such a poor approximation at 
long times. 3. The STA approximation for docs 
not satisfy either set of known sufficient conditions for 
symmetry exactly. However, in the limit in which the 
repulsive forces become hard sphere forces, it does sat- 
isfy one set. Thus for very strong repulsive forces it is at 
least numerically close to satisfying the symmetry condi- 
tions. 4. The RPMF approximation for satisfies the 
symmetry conditions exactly. 5. The hybrid approxima- 
tion discussed above does not satisfy either set of known 
sufficient conditions for symmetry exactly. 6. If a binary 
collision approximation satisfies either of the two known 
sets of sufficient conditions, its K approximation satis- 
fies the symmetry conditions exactly. Thus the K-MGBE 
and K-RPMF satisfy the symmetry properties exactly. 



APPENDIX B: TWO FORMS OF THE LEFT 
VERTEX IN THE MEMORY FUNCTION 

a. The left vertex in the MGBE and BCA. The rea- 
soning that leads to the derivation of the MGBE approx- 
imation in the context of the Mazenko theory and the 
BCA in the context of the graphical theory are very sim- 
ilar, but the results are slightly different. The nature of 
the difference is discussed in this appendix. 

The graphical theory contains a vertex Qi2(l; 1'2') de- 
fined as 

Oi 2 (l;i'2') 

= i J dl"d2"W 12 (l;l"2")K 2 (l"2"-l'2') (Bl) 
= (1 +1(1' 2')) V R u(12') • Vpd(ll') (B2) 

(See eq. (3) and Sec. Ill of Ranganathan and Andersen 42 
and Appendix A. 4 of Andersen4£ Here the 1 operator 
interchanges the two arguments in the subsequent ex- 
pression.) The first equality is a definition. The second 
equality represents an exact evaluation. Although this 
function has three arguments, they represent only two 
distinct particles, since the Dirac delta function in each 
term clearly indicates that the left argument corresponds 
to the same particle as one of the right arguments. So 
this vertex appears appropriate for a theory that includes 
binary collisions between particles. 

An exact evaluation of Wi2(l;l'2') gives two terms. 



The first term is 

W 12 (l; l^Ofirstterm = -(1 + X(l"2"))n(l>(2") 

x M m (l")M m (2")g(l"2")V P 6(l"l) . V rVmf (i"2") 

It clearly represents two distinct particles. The second 
term (not given here), which is very complicated, involves 
static correlations functions for three distinct particles. 

An exact evaluation of 14 7 i2(l;l'2') gives two terms. 
The first term is 

^i2(l;l'2') firs tterm = -(I + T(l"2"))n(l")n(2") 
x M m (l")M m 2")g(l"2")Vp<5(l"l) • V R v MF (l"2") 

This clearly represents two distinct particles. The second 
term (not given here), which is very complicated, involves 
static correlations functions for three distinct particles. 

An exact evaluation of Kiil" 2"; 1'2') also gives two 
terms. The first term is 

K 2 (l"2"; l'2') first term 
= (1+1(1"2"))S(1"1')5(2"2') 

n(l)M m (l)n(2)M m (2)g(12) 1 ' 

This term describes two distinct particles. The second 
term involves static correlations functions for more than 
two distinct particles. 

In the spirit of a binary collision approximation, it 
would be reasonable to approximate both W\ 2 and K 2 
by their first terms. If this is done and the results are 
used in (|B1[) . we get the following approximate result. 

Qi2(l;l'2') (B4) 
« (1+X(1 / 2 / ))Vrwmf(12')-Vp<5(110 

which is of the same form as the exact result but with 
the potential replaced by the potential of mean force. 

Mazenko's theory contains quantities analogous to W\ 2 
and K 2 . In the derivation of the MGBE, the contribu- 
tions that correspond to more than two distinct particles 
are not included in each of these two functions. This 
would be analogous to using eq. (|B4j) rather than (|B2jl . 
This is the origin of the fact that the final expression for 
the MGBE differs from the BCA in the way shown in 
eqs. and (|10[) . It is only in the factor of Qi 2 that 
restricting attention to two particle effects leads to dif- 
ferent results depending on whether the Q vertices or the 
separate W and K functions are regarded as the quan- 
tities that are to be expressed in terms of two particle 
contributions. 

b. Two comments. 1. The use of the exact result 
for Qi2 in the BCA is the origin of the fact that the 
BCA does not have the symmetry properties discussed in 
Appendix [A] The matter of whether an approximation 
satisfies the symmetry properties is determined by the 
relationships among all the approximations made rather 
than the accuracy of individual approximations. 2. This 
difference between the left vertices of the two theories 
(with the bare potential being replaced by a factor of the 
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potential of mean force) appears similar to differences 
noted in approximate kinetic theories of hard sphere s 34 ' 61 
and of plasmas . 62 ' 63 ! 64 ' 65 The plasma case may be related 
to the present situation but is more complicated in that 
both the left and right vertices are different in various 
versions of the theory. The similarity to the hard sphere 
case is superficial, since the latter involves having dif- 
ferent numbers of factors of g(d), the pair correlation 
function of hard spheres at contact, in the memory func- 
tion. A factor of the pair correlation function cannot 
arise from replacing the bare potential by the potential 
of mean force. The differences in the hard sphere kinetic 
theories appear more closely related to approximating K 2 
by the first term in eq. (IB3j) and then approximating the 
latter by its 'disconnected' approximation 

1*2\± & i r Z J first term disconnected 

(1+J(1"2"))(5(1"1')(5(2"2') 
~ n(l)M m (l)n(2)M m (2) 

which is equivalent to multiplying the approximate K 2 
by a factor of the pair correlation function g. 

APPENDIX C: NUMERICAL CALCULATION OF 
THE MATRIX ELEMENTS OF THE MEMORY 
FUNCTION AND CORRELATION FUNCTIONS 

The matrix elements of the memory function can be ex- 
pressed as an ensemble average over a two particle prob- 
ability distribution function as follows 

M flff (k,t) = (W^(k,t; 12)) f ' dr4wr 2 g(r) 

Jo 

where r c is the cutoff distance associated with the poten- 
tial governing the dynamics, that is V22, and the quantity 
being averaged is 

W)w(k,t; 12) = VVq 21 (R) • VMPi) 

x W Ql2 (|R Qb |) • [e-^Vh p (P a ) - e -* k - R >VMP&)] 

The meaning of phase points a and b are described in 
Set. IIII A[ The first term in brackets gives the self part 
of the memory function, while the second term gives the 
distinct part. Both terms together constitute the total 
memory function. The two particle distribution function 
required for calculation of the average (Wp^Ck, t; 12)) 

Pl 2 l(12) cx pM m (P 1 )M m (P 2 )d(R 1 )g(R 2 )e(r c - |R 2 |) 

The average of W over this distribution is calculated by 
the two particle trajectory method of Ranganathan and 
Andersenri 3 - In this method, the coordinates and mo- 
menta called 1, 2 are sampled from Pl 2 l and used as ini- 
tial conditions for a two particle trajectory calculation 
from which the W quantities are calculated as functions 
of time. For most of the approximations, the trajecto- 
ries are calculated using Hamiltonian mechanics of the 
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FIG. 10: A schematic representation of the averaging proce- 
dure for the matrix elements of the memory function. 



two particles interacting via the potential V-22- For the 
K approximations, the dynamics are supplemented by 
an algorithm that includes the effect of the Lbgk in 
cq. (fTTj) . The algorithm gives each particle, at randomly 
chosen times, a new momentum chosen from a Boltz- 
mann distribution, without changing its position. The 
time at which these interventions occur for each particle 
are uncorrelated with each other and occur with an av- 
erage frequency v for each particle, where v is given by 
eq. ([TO]). 

The matrix elements of the memory function can be 
divided into four components: self real, distinct real, self 
imaginary, and distinct imaginary parts. Each part is 
computed for each matrix element in a calculation that 
takes advantage of a number of properties derived from 
the Hermite polynomial basis functions, the two parti- 
cle problem, and the known t = values of the matrix 
elements themselves. 

The basis function h^(P) is a constant. Since the gradi- 
ents of the Hermite polynomial functions with respect to 
momenta, not the Hermite functions themselves, appear 
in the expression for the matrix elements of the memory 
function, Mg- for all V and M-g for all ft are identically 
zero for all time and need not be computed numerically. 

The two particle trajectory method can be mapped 
onto a highly symmetric two particle scattering prob- 
lem. The symmetry relations of most concern are the 
following. 1. If fi z + v z is even, then the self and distinct 
imaginary parts are identically zero for all times. 2. If 
\i z + v z is odd, then the self and distinct real parts are 
identically zero for all times. See Ref. |48| for the proofs 
of these properties. 

Figure [TO] is a schematic representation of the averag- 
ing procedure used to calculate the matrix elements of 
the memory function. For each pair of indices /2z7, the 
following calculation scheme was used. 1. A set of initial 
conditions is randomly sampled from the two particle 
probability distribution and a molecular dynamics 
trajectory is computed by solving Hamilton's equations 
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of motion. 2. This single trajectory is used to generate 
fifteen symmetry related two particle trajectories based 
on the Dih symmetry of the two particle system. The 
values of the various parts of Wp,p are calculated for each 
of the sixteen trajectories. This corresponds to gener- 
ating numbers that appear in the rightmost column in 
Fig. [TUJ The sixteen sets of numbers are averaged to- 
gether to give the first type of average for the various 
parts of Wait- The values are entered into the middle 
column of the figure. 3. Steps 1 and 2 are repeated for a 
total of N T times. The N T results are averaged to give 
a second type of average, which is entered in the left col- 
umn of the figure. 4. Steps 1-3 are repeated Ne times 
to yield a set of Ne values of the second type of aver- 
age. These Ne values are averaged together to give a 
third type of average, however they are not given equal 
weight as in steps 2 and 3. In this case, we use a con- 
strained reweighting method, based on a maximum likeli- 
hood principle, that adjusts the weights assigned to each 
of the Ne values such that the correct results are ob- 
tained for the calculated averages of the matrix elements 
at t — 0. The result of this process is a set of values for 
the real and imaginary parts of the self and distinct parts 



of (WfCff). The values automatically satisfy the symme- 
try properties mentioned above, because of the way the 
symmetry related trajectories were used, and they auto- 
matically have the correct values at t = 0, because of the 
reweighting method. The latter automatically reduces 
the statistical error for small nonzero times as well. De- 
tails of this method can be found in Ref. |48|. The self 
parts of the average are used in eq. (fD2")l to give the de- 

A (s) 

sired matrix elements M~J. 5. The correlation functions 

of interest in this work, Ci s ] (k, t) and C^) (k, t) , are com- 
puted by numerically solving the kinetic equation using 
the Euler method with the appropriate matrix elements 
of the memory function. 

Steps 1-5 are repeated N s times to generate N 8 statis- 
tically independent results for the functions Mn^it) and 

for CL s J(t) and Ci s ; (t). These N s results for each func- 

00 w kk w _ s 

tion are averaged, and statistical error bars are computed 
in the usual way. The results presented in this paper 
are these average quantities. The values of Nt, Ne and 
N s used in this work are N T = 100, N E = 1000, and 
N s = 10. 
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